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Abstract 

The present work deals with smooth-muscle-tissue material- 
model development and its use in finite-element 
computational analyses. In the first portion of the work, one 
of the smooth-muscle tissues reported in the open literature 
is critically assessed. Then, the model is calibrated using 
experimental data published in the open literature, for 
isometric (constant length) and isotonic (constant holding 
stress) uniaxial-tensile-tests of swine carotid artery. 
Predictions of this model under uniaxial loading conditions 
are then investigated and their physical soundness is 
established. Next, the model is implemented into a user- 
material subroutine which enables, under general 
deformation conditions, determination of the spatial 
distribution and temporal evolution of the material stress 
state and material state-variables. Lastly, the user-material 
subroutine is linked with a commercial finite-element 
program and used to analyze a simple smooth-muscle 
contractile problem involving: (a) a tubular structure 
surrounded with smooth-muscle tissue (with muscle-cell 
longitudinal directions aligned in the tube circumferential 
direction); and (b) an intra-luminal content in the form of a 
prolate spheroidal, the major axis of which coincides with 
the tube axis. Activation/contraction of the smooth-muscle 
tissue triggered by the tube/spheroidal contact stresses has 
been found to result in a forward momentum/motion of the 
intra-luminal content. 
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Introduction 

Within the present work, the contractile response of 
smooth muscle tissue is analyzed computationally and 
the response is correlated with the basic architecture/ 
microstructure of this tissue. Hence, the main aspects 
of the present work include: (a) basics of smooth 
muscle tissue; (b) prior computational and experimental 
work dealing with phenomena and processes 
accompanying smooth muscle contraction; and (c) 
multi-scale hierarchical structure of smooth muscle 


tissue. A brief overview of these aspects of the 
problem at hand is presented in the remainder of this 
section. 

77; c Basics of Smooth Muscle Tissue 

The defining features of smooth muscle tissue are: (a) 
unlike skeletal muscle tissue which can undergo only 
voluntary contractions, smooth muscle tissue 
undergoes involuntary contractions; (b) its contraction 
is associated with energy consumption which is ca. 
0.2-1% of that associated with skeletal muscle (Brophy, 
2000; Paul, 1990; Walker et al., 1994). As a result, 
smooth muscle tissue is capable of undergoing 
prolonged contractions without experiencing fatigue 
(in contrast to skeletal and cardiac muscle) (Morgan et 
al., 1989); (c) while smooth muscles are capable of 
delivering comparable isometric (i.e. constant muscle 
length) force per unit cross-sectional area as skeletal 
muscle, the rate of smooth muscle contraction is one to 
two orders of magnitude lower (Bitar, 2003; Walker et 
al., 1994); (d) smooth muscle tissue is typically found 
in organs which experience continuing contraction 
and relaxation such as the stomach, the intestines, the 
urinary bladder, the airways and blood vessels; and (e) 
the function of the smooth muscles is generally 
dependent on the system/organ within which they 
reside. For example, within the vascular system, blood 
flow is regulated by the smooth muscle tissue located 
within the arterial wall. On the other hand, within the 
digestive system, smooth muscles govern the mixing 
and transport of intraluminal contents and thus, 
ensuring: (i) efficient digestion of food; (ii) progressive 
absorption of nutrients; and (iii) evacuation of residues. 

Prior Studies of Smooth Muscle Contraction 

The physiology of smooth muscle contraction is highly 
complex, and involves interaction among a number of 
(not yet well-understood) electrical, biochemical and 
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mechanical phenomena and processes. Since 
abnormalities in smooth muscle contraction can have 
severe consequences and may cause diseases such as 
hypertension and asthma, smooth muscle contractions 
have been the subject of many experimental and 
theoretical studies. 

Among the experimental studies dealing with the 
smooth-muscle contraction, the following appear the 
most noteworthy and relevant to the present work: 
smooth muscle from (a) guinea pig taenia coli (Arner, 
1982; Lofgren et al., 2001; Peterson, 1982); (b) porcine 
carotid artery (Hai and Murphy, 1989; Kamm et al., 
1989; Rembold and Murphy, 1988, 1990; Roy et al., 
2005; Silver et al., 2003; Singer and Murphy, 1987); (c) 
porcine coronary artery (Makujina et al., 1995); (d) 
porcine tracheal (Herrera et al., 2005); (e) canine 
carotid artery (Takamizawa and Hayashi, 1987); (f) 
ferret aorta (Jiang and Morgan, 1989); (g) bovine 
tracheal (Tang et al., 1992); (h) rat aorta (Tosun et al., 
1997); and (i) rat pulmonary artery (McIntyre et al., 
1996). 

As far as the theoretical studies are concerned, they 
could be broadly divided into two groups: (a) those 
dealing with the mechanism of smooth muscle 
contraction (Fay and Delise, 1973; Herrera et al., 2005; 
Miftakhov and Abdusheva, 1996; Rachev and Hayashi, 
1999; Rosenbluth, 1965; Yang et al., 2003a,b); and (b) 
those focusing on the development of the appropriate 
continuum-level material constitutive model (Gestrelius 
and Borgstrom, 1986; Lee and Schmid-Schonbein, 
1996a,b; Stalhand et al., 2008; Zulliger et al., 2004). 

The work reported in the present manuscript is most 
closely related to the work of Kroon (2010) who 
proposed a theoretical model for the constitutive 
behavior of smooth muscle tissue, undergoing large 
deformations. The model includes both the contributions 
of the calcium-activated contraction and the passive 
deformation-resistance to the constitutive response of 
smooth muscle tissue. The active portion of the 
material response is modeled using the concepts of: (a) 
contractile units consisting of thick myosin filaments 
and thin actin filaments; and (b) filament sliding. The 
passive portion of the material response is of a 
viscoelastic character and includes the contribution of 
the remaining intracellular content as well as of the 
extracellular matrix. 

Multi-scale Hierarchical Structure of Smooth Muscle 
Tissue 

As schematically shown in Figure 1(a), smooth 
muscles possess a layer structure with each layer/sheet 


containing spindle-shaped cells (each having a single 
nucleus). Typical dimensions of a smooth muscle cell 
are: length = 200-300 //m and width = 5-15 /an (Bitar, 
2003). 



FIGURE 1(a) A SCHEMATIC OF THE SMOOTH-MUSCLE-TISSUE 
SHOWING LABELED KEY COMPONENTS; AND (b) A 
SCHEMATIC OF THE CLOSE-UP OF A CONTRACTILE 
APPARATUS UNIT 

Based on their structural/mechanical function, the 
following main components of a smooth muscle cell 
are typically identified: 

(a) the "contractile apparatus," which consists chiefly 
of thick myosin filaments and thin actin filaments. It 
should be noted that in vivo, these filaments are quite 
labile, and may undergo major structural changes/ 
reorganizations. Typically, the thin filaments are more 
abundant, with the ratio of the thin to thick filament 
numbers as high as 15 in the case of vascular smooth 
muscle. In general, the myosin/actin contractile- 
apparatus assemblies are found to be approximately 
parallel to the cell longitudinal axis (Bitar, 2003; 
Herrera et al., 2005; Hodgkinson et al., 1995; Kuo and 
Seow, 2004; Seow and Par, 2007). Often bundles of 
closely-spaced parallel thick myosin filaments are 
found to terminate at the same longitudinal location of 
the cell. Thick filaments in vascular smooth muscles 
with a typical length of 2.2 jum are longer than their 
counterparts in skeletal muscle. It should be further 
noted that to achieve a large contraction stroke and its 
rate, the length of the contractile unit is quite short in 
comparison to the cell longitudinal direction (i.e. the 
number of contractile units across the longitudinal 
dimension of the cell is quite large). 

The architecture of the smooth muscle cells is 
consistent with the so-called "side-polar conformation" 
of the thick myosin filaments (the conformation within 
which cross-bridges on the opposite sides of the 
filament have opposite orientations), and with the 
actomyosin side-polar interaction-based filament 
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sliding model of smooth-muscle contraction (Herrera 
et al., 2005; Hodgkinson et al., 1995). The sliding 
filament theory (Huxley, 1953, 1957; Huxley and 
Niedergerke, 1954; Huxley and Hanson, 1954) 
postulates that the force required for muscle contraction 
is produced by advancement of the myosin 
heads/cross-bridges along the thin actin filaments. 
Figure 1(b). As seen in Figure 1(b), the side polar 
model of actomyosin interaction results in the actin 
filaments that interact with different ends of the thick 
myosin filament, sliding in opposite direchons. 

(b) the "cytoskeleton," which contains actin filaments 
for passive structural support of the cell, as well as 
intermediate/connecting filaments; 

(c) "dense bodies," which contain a-actinin (an actin 
cross-linking protein) and function as the focal contact 
points between the contractile apparatus and the 
cytoskeletal fibers within the cytoplasm. In addition to 
being a focal point to the actin thin filaments 
(Hodgkinson et al., 1995), dense bodies are sometimes 
seen to also act as an anchor to the intermediate 
filaments. Generally, about 80% of the cell interior is 
occupied by dense bodies and contractile filaments 
(Bitar, 2003); and 

(d) "dense plaques," the inter-cell focal contact points, 
located throughout the plasma membrane, which also 
act as anchors to the actin thin filaments and 
cytoskeletal intermediate filaments. Since contact 
interactions between pairs of opposed adjacent dense 
plaques located on neighboring cells enable force 
transmission across cell boundaries, the contractile 
apparatus effectively crosses cell boundaries. 

Smooth muscle tissue and their cells are typically 
analyzed computationally within a continuum 
framework. That is, structural details of the tissue/cells 
are smeared out into a homogeneous (often isotropic) 
continuum. The smooth-muscle-tissue model used and 
extended in the present work falls into this category of 
material models. However, there are finer-scale non- 
continuum material models for the components of the 
smooth muscle cells. Since, as established above, about 
80% of the cell interior is occupied by contractile 
filaments (including dense bodies) (Bitar, 2003), and 
these finer-scale material models typically focus on 
capturing various aspects of constitutive response of 
the contractile apparatus and its components. As seen 
in Figure 2, the following additional length-scales/ 
constitutive models are typically used to capture the 
behavior of contractile elements and their components: 
(a) a structural element which models the contracting 


apparatus as a structural unit consisting of truss-type 
and beam-type elastic and dissipative elements 
(Grujicic et al., 2008ab, 2009ab, 2012); (b) a meso- 
length-scale model within which beads (hypothetical 
particles composed of a large number of atoms/ions) 
and coarse-grained molecular statics/dynamics are 
used to model details of the myosin-head/actin- 
filament interactions (Chu and Voth, 2005, 2006); and 
(c) atomic-scale modeling of various aspects of 
microstructure and properties of different components 
of the contractile apparatus, e.g. G-actin monomers, F- 
actin polymers/filaments, etc. (Chu and Voth, 2005; 
Pfaendtner et al., 2010). As mentioned above, the 
present work utilizes the concept of a continuum 
material to treat the contractile and the accompanying 
passive response of the smooth muscle tissue and cells. 


Continuum/Tissue-Scale Contractile-Unit Scale 



\/ 


Molecular-Scale Filament Meso-Scale 



FIGURE 2 MULTI-LENGTH-SCALE HIERARCHICAL 
STRUCTURE OF THE SMOOTH-MUSCLE-TISSUE 


Main Objective 

The main objective of the present work is to extend the 
smooth-muscle-tissue material constitutive model 
proposed by Kroon (2010) so that the model can be 
applied within a finite-element computational 
framework. This extension then enables modeling of 
the contractile and passive responses of a simple organ, 
such as a short section of the (small) intestine. As will 
be shown later, to obtain the extended version of the 
model proposed by Kroon, a material user subroutine 
had to be developed and linked with the general- 
purpose finite element code. 

Paper Organization 

An overview of the basic biomechanics of smooth 
muscle contraction is provided in Section II. Chemistry 
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aspects of the smooth-muscle model proposed by 
Kroon (2010) are discussed in Section III, while the 
mechanical aspects of the same model are assessed 
critically in Section IV. Application of the Kroon 
model to the cases of isometric (constant muscle length) 
and isotonic (constant muscle force) uniaxial tensile 
tests is presented and discussed in Section V. The 
procedure used to parameterize the model for the case 
of swine carotid artery smooth-muscle-tissue is 
described and the resulting model parameters listed in 
Section VI. Validation of the model as well as an 
analysis of its key predictions are presented in Section 
VII. A detailed description of the extension of the 
material model to the finite-element computational 
framework and its use to mimic the response of a 
simplified organ (a short section of the small intestine) 
is provided in Section VIII. Key conclusions resulting 
from the present work are summarized in Section IX. 

Biomechanics of Smooth Muscle Contraction 

Structural and Mechanical Aspects 

Since the basic structural and mechanical aspects of 
smooth-muscle tissue were reviewed in the previous 
section, they will not be repeated here. Instead, the 
reader is directed towards reviewing Figures l(a)-(b) 
for appropriate details. 

Biochemical Aspects 

While either mechanical, electrical or chemical stimuli 
may trigger smooth muscle contraction, in each case 
muscle contraction is preceded by an increase in the 
concentration of intracellular Ca 2+ ions, denoted by 
[Ca 2+ ] (from ~0.1 pM, corresponding to the relaxed 
state of the cell, to ~1 ,uM in the contracted state of the 
cell). 

An increase in the intracellular [Ca 2+ ] sets off the 
following sequence of biochemical reactions: (i) Ca 2+ 
binding to calmodulin is first triggered; (ii) myosin 
light chain kinase (MLCK) is then activated by the 
calcium/calmodulin complex which binds to it (Jaggar 
et al., 2000); and (iii) activation of MLCK in turn 
activates (independently) each myosin head. Activation 
of the myosin heads enables cycling/advancement of 
the cross-bridges (connecting the thick myosin 
filament to the thin actin filament) along the actin 
filaments and, in turn, generation of the force required 
for muscle-cell contraction/ shortening. 

The cycling of cross-bridges itself involves a sequence 
of coupled/correlated biochemical and mechanical 
phenomena and processes, including: (i) an adenosine 


triphosphate (ATP) molecule (the source of energy for 
the contraction) binds to the myosin head, causing its 
phosphorylation; (ii) the new conformation of the 
myosin head results in its rapid detachment from actin 
(if the two were bound); (iii) next, hydrolysis of ATP 
takes place, which produces adenosine diphosphate 
and inorganic phosphate (Pi); (iv) a significant change 
in the myosin-head conformation is associated with 
the ATP hydrolysis; (v) the myosin-head in the new 
conformation attaches, through the formation of cross- 
bridges, to the adjacent actin filament; (vi) next, the Pi 
dissociates from the myosin head, allowing the release 
of energy associated with ATP hydrolysis. This, in 
turn, leads to cycling/peddling of the cross-bridges 
and to relative sliding between thick myosin filaments 
and thin actin filaments. To prevent back-sliding of the 
actin filaments, cross-bridges do not cycle in unison. 
Examination of the sequence of biochemical and 
mechanical processes just described reveals that: (a) 
cell contraction is dependent on the hydrolysis of ATP, 
the process which provides the necessary chemical 
energy and its conversion into mechanical work; and 
(b) the force generated during muscle contraction is a 
function of the number of cross-bridges connecting the 
thick myosin filament to, and advancing along, the 
thin actin filaments (Bitar, 2003). 

Under sustained (isometric) contraction conditions, 
the mechanism of cross-bridge cycling is believed to 
be different from the one described above. In this case, 
initial muscle contraction is still activated by an 
increase in the cytosolic [Ca 2+ ]. However, once the 
maximum contraction of the muscle is attained, the 
cytosolic [Ca 2+ ] decreases (without an accompanying 
decrease in contractile stress). In this state, the cross- 
bridges become dephosphorylated while remaining 
attached to the actin filaments. This state is commonly 
referred to as the "latch" state (Dillon et al., 1981) since 
the cross-bridges either cycle infrequently or do not 
cycle (Bitar, 2003). In the latch state, cross-bridges 
remain under tension while maintaining force in 
sustained muscle contraction. This is accomplished, 
however, at a level of ATP consumption almost an 
order of magnitude smaller than the level that would 
be required under normal cross-bridge cycling 
conditions (Clark and Pyne-Geithman, 2005). 

Chemical Model of Smooth Muscle Contraction 

As mentioned earlier, according to the sliding filament 
theory (Huxley, 1953, 1957; Huxley and Niedergerke, 
1954; Huxley and Hanson, 1954), muscle contraction is 
a result of the relative sliding of myosin and actin 
filaments within the contractile apparatus. The cross- 
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bridges play a key role in the biochemistry and 
mechanics of smooth muscle contraction, as mentioned 
earlier. The myosin/actin contractile apparatus is 
simulated in the present work using the model of Hai 
and Murphy (1988), which postulates the existence of 
four states (named si, s^, S3, S4 respectively) of the cross- 
bridges. The existence of the four states is a result 
of the actions of MLCK and myosin light chain 
phosphatase on both free and actin-attached cross- 
bridges. According to this model, only phosphorylated 
cross-bridges are capable of attaching to actin 
filaments. The four states can be briefly described as 
follows: (a) when the cross-bridge is in state si, it is 
dephosphorylated and passive and, hence, unable to 
attach to the actin filament; (b) in state si, the myosin 
head is phosphorylated by the action of the MLCK, i.e. 
an ATP molecule is attached to the cross-bridge. 
However, while the myosin head in the S 2 state is 
capable of attaching to the actin filament, this 
attachment has not occurred yet; (c) upon attaching to 
the actin filament, the cross-bridge assumes state S3. It 
is the cycling of cross-bridges between states S 2 and S3 
which creates the muscle contractile force; and (d) 
while remaining attached to the actin filament, some 
myosin may become dephosphorylated. This (S4) state 
of the cross-bridge, as mentioned earlier, is commonly 
referred to as the latch state. While myosin heads in 
state 4 do not cycle any further and, hence, do not 
contribute to relative filament sliding, they generate a 
passive resistance force within the contractile 
apparatus (and, thus, help slow down the rate of 
detachment of the filament heads from the actin 
filaments). Based on the description of the four 
myosin-head states, it is clear that only states S3 and S4 
can generate/transmit force between the myosin and 
actin filaments. 


The contractile-apparatus four-state model of Hai and 
Murphy (1988) includes the following basic assumptions 
and simplifications: (a) the total number of myosin 
heads that may potentially bind to an actin filament, 
mtot, is assumed to be constant; (b) the fraction of 
these myosin heads residing, at any instant of time, in 
state Si (z=l-4) is denoted as m and the four n,'s act 
as smooth-muscle-contraction chemical-model state 
variables; (c) evolution of the chemical-model state 
variables is assumed to be given by a first-order 
chemical-reaction model in the form 


dn 

dt 



k 2 0 

-k 2 -k 3 k 4 

k 3 k 4 k 5 

0 k 5 


0 

k„ 

-k 6 - k 7 


n = Kn (1) 


where n — \n x n 2 n 3 n 4 ] (the superscript T denotes a 
transpose operator) is a state-variable column vector 
while ki through h denote rate constants associated 
with the transitions between different cross-bridge 
states (K is the rate-constant matrix); (d) Eq. (1) is 
subjected to the conservation law which stipulates that 
the sum of the four n’s, at any time instant, is equal to 
1. It should be noted that, as mentioned earlier, muscle 
contraction is triggered by an increase in the cytosolic 
[Ca 2+ ]. Since this increase requires an increase in the 
value of m, it causes an increase in the rate constants k\ 
(controlling the si- S 2 transition) and kt, (controlling the 
S4 — S3 transition). To help clarify the transition between 
the four cross-bridge states, a schematic of these states 
and their inter-relationship is provided in Figure 3. 



FIGURE 3 FOUR STATES OF THE CROSS-BRIDGES PROPOSED 
BY HAI AND MURPHY (1988), AND LABELS FOR THE RATE 
CONSTANTS RELATED TO THE TRANSITIONS BETWEEN 
DIFFERENT STATES 

Mechanical Model of Smooth Muscle Contrac- 
tion 



FIGURE 4 RHEOLOGICAL MODEL OF THE SMOOTH-MUSCLE 
TISSUES PROPOSED BY KROON (2010). PLEASE SEE TEXT FOR 
DETAILS 
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According to the preceding discussion and Figure 1(a), 
contractile-apparatus filaments are approximately 
aligned in the cell longitudinal direction and transcend 
the cell boundaries. In the lateral direction, these 
filaments are closely-spaced and parallel, with the 
myosin heads located at fairly regular longitudinal 
intervals. Taking this consideration into account. 
Kroon (2010) proposed a muscle-cell mechanical 
model which contains the following elements/ 
assumptions: (a) contractile-apparatus filaments are 
assumed to all be parallel, aligned with the cell 
longitudinal direction and surrounded by a passive 
matrix. The latter contains: (i) the remainder of the cell 
interior (composed of cell membrane, cytosol, passive 
networks of actin, intermediate filaments, cell nuclei); 
and (ii) extracellular matter typically composed of 
elastin and collagen; (b) to capture this behavior, the 
corresponding rheological model is constructed which 
contains three parallel branches, one " active " 
(generates the contractile force) and two passive 
(Figure 4), aligned in the muscle-cell longitudinal 
direction; (c) the constitutive response of the 
rheological model displayed in Figure 4 is described 
using the large-deformation formalism and the 
concept of the strain-energy (density) function W as: 

^ = ^act+^pass (2) 

where subscripts act and pass are used to denote the 
active and passive components of W; (d) to account for 
the potentially large deformation and rotations 
accompanying muscle contraction, a single Cartesian 
coordinate system defined in terms of three basis 
vectors e 2 , e 3 is first introduced. Then, the position 


vectors of a material point in the reference 
configuration Do and the current configuration D are 
defined as X = A/e, and x = x i e i , respectively; (e) the 
kinematics of deformation is next described using the 
deformation gradient F : 

F = dx/dX, (3) 

right Cauchy-Green deformation tensor, C : 


C = F r F 


(4) 


(5) 


and the Green-Lagrange strain E : 

E = — (C - 1) 

2 ' 

where I is the second-order identity tensor; (f) the 
stress state is defined by the second Piola-Kirchhoff 
stress S as: 

s _^_ 2 ^_ 2 ^ +2 y<w 

5E DC DC dC p ' l j 


by the first Piola-Kirchhoff stress P as: 

p = fs = f(s, i ,+s [ ,J=p„+p i 


pass 


(7) 


and by the Cauchy (true) stress a as: 

<T = FSF r /detF ( 8 ) 

where det denotes a determinant operator. 


Active Part of the Mechanical Response 


As shown in Figure 4, the contractile apparatus is 
modeled using one of the three branches, i.e. the active 
branch. The main characteristics of the active branch 
can be summarized as follows: (a) this branch is 
assumed to consist of an ‘‘active dashpot" (a contractile- 
force-generating time-dependent element) and a 
passive spring. The active dashpot is used to describe 
the cycling behavior of the cross-bridges which results 
in relative filament sliding, while the passive spring is 
used to describe the elastic resistive response of the 
cross-bridges. It should be noted that the myosin thick 
filaments and actin thin filaments, in the present 
contractile-apparatus formulation, are treated as rigid 
bodies. In other words, the entire filament-sliding and 
passive-resistance response of the contractile element 
is assumed to be localized within the myosin heads; (b) 
the active branch is aligned with direction M (defined 
as a unit column vector) in the reference configuration 
Do; (c) stretching of the (entire) structural element in 

the filament direction (i.e. direction M ), 'V , is defined 

as %f = M CM . The effect of filament contraction is 

depicted in Figure 5, which shows that the distance 
between two adjacent dense bodies (bridged by a 
contractile-apparatus) has been reduced from its initial 

value L cu (in the reference configuration) to ^f^cu (in 
the current configuration); (d) within the active branch, 

7/ is decomposed multiplicatively as = A.Vty f 


where ^fc accounts for the relative sliding between 
myosin and actin filaments ( fc stands for filament 
contraction) and A, e accounts for the elastic deformation 
of the cross-bridges; (e) the following ‘‘active strain 
energy function" is utilized (Fung, 1970): 


!^act = ^p(«3 + **)(% - 1 ) 2 


Pf / \ 

= ^-(«3 +l h) 


f X 2 

MCM 


(9) 


where: (i) Mf is a stiffness parameter which is 
proportional to the stiffness of cross-bridges and to the 
number of parallel contractile filaments per unit area 
of the cell transverse section; (ii) the fact that only 
cross-bridges in states S3 and S4 contribute to a force 
generation/transmittal within the active branch is 
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reflected by the use of the (773 + / 7 4 ) term; (iii) cross- 
bridges in states S 3 and S 4 are assumed to possess the 
same elastic stiffness; (iv) the total number of myosin 
heads, at a given level of [Ca 2+ ], participating in the 
construction of a contractile apparatus, niiot , remains 
constant; and (v) the functional form of '/act ensures 
that, in the relaxed state, not only the strain energy but 
also the three aforementioned stress measures are zero; 
and (f) examination of the active-branch mechanical 
model shows that it is associated with four state 
variables (i.e. variables which define the current state 

of the contractile-apparatus): « 3 , n 4 , and ^fc . To 
compute the instantaneous values of these state 
variables, the corresponding evolution equations must 
be defined and integrated in time. 




FIGURE 5 A SCHEMATIC OF THE FILAMENT SLIDING 
MECHANISM FOR SMOOTH MUSCLE CONTRACTION 
SHOWING: (A) THE RELAXED; AND (B) THE CONTRACTED 
CONFIGURATIONS OF A SINGLE CONTRACTILE- APPARATUS 

1) Evolution Equations for « 3 and n 4 

Evolution of « 3 and n 4 has been previously 
defined using Eq. (1). 

2) Evolution Equation for 'y 

2 T 

Since Xj- =M CM, its evolution is defined by the 
temporal evolution of C (i.e. of F )• 


3) Evolution Equation for X fc 

Kroon (2010) proposed the following evolution 

equation for ^ fc : 


n^fc ~ p C c 


dX 


'/c 


( 10 ) 


where P cc denotes the active stress generated by 
the active dashpot in response to an increase in 
[Ca 2+ ], q is a viscous-damping coefficient associated 
with the cross-bridge-cycling/fiber-contraction 
process, and the term j s the elasticity- 

type myosin-head resistance to fiber contraction. 
To construct the appropriate expression(s) for P cc , 
Kroon made the following observations: (a) since 
P cc is generated as a result of muscle-cell 
contraction, this term is of a compressive-negative 
character; (b) only phosphorylated and attached 
cross-bridges in state S 3 contribute to filament 
contraction; (c) when P cc does not possess a 
sufficiently high magnitude to cause further 
muscle contraction, but is comparable to 
S V^act / 5A /C ; dephosphorylated attached cross- 
bridges in state S 4 can help prevent backsliding in 
the contractile filaments; (d) when the elastic- 
resistance term -di// act /dXj- c exceeds the magnitude 


of P cc , the dephosphorylated cross-bridges become 
detached from the actin filaments and no longer 
contribute to the load transfer within the active 
branch. While cross-bridges in state S 3 may also get 
detached, their phosphorylated character allows 
them to quickly reattach to the contracting filament 
and to slow down the rate of filament 
backsliding/extension. To capture this behavior of 
the active branch. Kroon proposed the following 
expressions for P cc : (i) for the contracting state of 

the filament ( X/ c < 0 ), within which the magnitude 


of the active stress exceeds the elastic resistance: 

C^act 


P cc=~ K 3 « 3 > for k " 3 « 3 > - ■ 


dX 


■fc 


( 11 ) 


(ii) for the intermediate state of the filament 
( Xf c - 0 ), which is characterized by nearly zero 


changes in V 

p g^act 

r CC ~ 0 5 

sx fc 


for /c 3 « 3 < 


d^act 

dX fc 


< K^n 3 + at 4 «4 


(iii) for the backsliding/extension state of the 
filament (X/ c > 0 ): 


-/c 3 « 3 , for /c 3 « 3 + K 4 n 4 < 


d^act 

5X fc 


(13) 


where /r 3 and /r 4 reflect the intensity/strength of 
the filament-contraction driving force and the 
filament-contraction elastic resistance, respectively. 


Specification of l^act via Eq. (9) and the functional 
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relationship given in Eq. (6) enables determination 
of the active second Piola-Kirchhoff stress as: 


S =2 


h f (n 3 +n 4 ) 


dC 




M T CM 


-1 


M ® M (14) 


where ® denotes the dyadic product of two vectors. 
4) Summary of the Model Parameters 


Examination of Eq. (14) shows that the active part 
of the smooth-muscle-tissue mechanical model 

contains four material parameters: , K 3 , K 4 and 

ri- 


Passive Part of the Mechanical Response 

1) Decomposition of the Deformation Gradient 

The passive part of the smooth-muscle material is 
treated as being isotropic linear-viscoelastic and is 
modeled using a generalized Maxwell model 
consisting of two parallel branches. Figure 4. In 
addition, in accordance with the general practice, it 
is assumed that only the isochoric (i.e. volume- 
conserving) part of deformation can be viscous (i.e. 
time-dependent). In other words, the volumetric 
part of deformation is assumed to be purely elastic. 
Consequently, the total deformation gradient is 
decomposed multiplicatively as: 

F = / 1/3 F (15) 

where J = detF is a measure of the deformation- 
induced volume change, and F defines the 
isochoric part of the deformation gradient (i.e. 
detF = 1 )• 

Since stresses arise solely from the elastic part of 
the deformation, the isochoric part of the 
deformation must be further decomposed into its 
viscous and elastic parts as: 

F = F e F v (16) 

where subscripts e and v are used to denote, 
respectively, the elastic and the viscous parts of a 

quantity, and both F e and F v are isochoric. 

Using the relation given by Eq. (4), the total 
isochoric, elastic and viscous types of the right 
Cauchy-Green deformation tensors can be defined 
as: 

C = F r F (17) 

C.=F/F. (18) 

C v = f/f v (19) 


defined as the sum of three (i.e. volumetric, elastic- 
isochoric and viscous) terms as: 

^pass =iPpa»(^C,C T ) 

( 20 ) 

= U(J) + y„( C) + ¥v (C,C y ) 

where the three parts (It, ^ c0 , 1/7 '' ) of the passive- 
strain-energy function are, respectively, defined as: 



where K, the bulk modulus, and - Ue , the steady- 
state "shear-spring constant " of the passive part of 
the muscle, are associated with the bottom (single- 

spring branch) in Figure 4, and C, = J~ is the third 

invariant of C . Pt l e denotes the shear-spring 
constant within the second/viscous branch while • 
represents the inner/scalar product of two second- 
order tensors. 


Examination of Eqs. (21)-(23) reveals that the 
deformation state of the passive part of the smooth- 
muscle-tissue mechanical model is defined by C 
and C v (and their invariants). 

3) Passive-Stress Components 

Via the use of Eqs. (6) and (21)-(23), the 
corresponding components of the passive second 
Piola-Kirchhoff stress can be derived as: 

S vol =K(J-l)JC~ l (24) 

S 00 =// e / 3 - 1/3 ^I-i(C:I)C“ 1 j (25) 

S v = Pp e If m [c~ X ~(C : C; 1 )^ 1 j (26) 

while the total passive second Piola-Kirchhoff 
stress is given as: 

O _ ry ^ W paSS 

pass - dc 

= 2 dUfJ) | 2 dyr x | 2 dy/ y (27) 

SC dC dC 

= S vol +S 0 O +S v 

The corresponding expressions for the first Piola- 
Kirchhoff stress and for the Cauchy stress can be 
obtained using Eqs. (7)-( 8), respectively. 

4) Evolution Equation for S v : 


2) Passive-Strain-Energy Function 

The passive part of the strain-energy function is 


Examination of Eq. (23) shows that as C v relaxes 
towards C = / 3 C , Vv , and hence, S v vanish. In 
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order to capture time evolution of the viscous 
second Piola-Kirchhoff stress and to obtain a 
reasonable agreement with the experimental data. 
Kroon (2010) proposed the following functional 
relationship: 


f dS v 

\ 


( 3S V , ^ 

k dt 

C=const / 

V 

Uc/ V J 


= -I(5v > .. + « 2 sgn(5 v ..)S v 2 .. + « 3 5 l 3 J 

where r is the relaxation time. Typically, only the 
first term associated with the exponential 
relaxation of the form exp(-f/r) r where t is time, is 
used. However, in order to obtain a good 
agreement with the experimental results. Kroon 
found it necessary to introduce the two additional 
non-linear terms in Eq. (28) with u 2 and a 2 being the 
corresponding coefficients. It should be noted Eq. 
(28) implicitly defines C v and, hence, the evolution 
of C v . 

5) Summary of the Model Parameters 

Examination of Eqs. (15)-(28) shows that the 
passive part of the smooth-muscle model includes 
the following model parameters: K , n e , P, T - a i 
and a 2 . 


Thermodynamic Validity of the Model 

Kroon (2010) also established that the model is 
thermodynamically valid, i.e. the internal energy is 
always non-negative, provided the following 
two conditions are satisfied: (a) y>0 ■ and (b) 

(-a^/ac v ):C v > o. 

Model Reduction for Uniaxial Tensile Loading 
Case 

Since the muscle-cell material-model parameters, 
summarized in the previous sections, are typically 
determined using laboratory uniaxial tensile tests, it is 
beneficial to recast the model into its uniaxial tensile 
form. This is done in the present section. It should be 
noted that, for convenience, the filaments are aligned 
in the e 3 direction, i.e. M = e 3 . 

1) Deformation Variables 

(a) Deformation gradient- F = diag{\,X} ,T 3 ) where 
diag denotes a diagonal second-order tensor, and 
4's(/ = l-3) are the three principal stretches; (b) 
Viscous deformation gradient- F v = diag(/L v j , X v j , X v ^ ) ; 
(c) Due to symmetry/trans verse-isotropy, \ - T 2 


and T v j = A v >2 ; (d) Due to the isochoric character of 
F v , ^v,l^v,2^v,3 ~ 1 or ^v,i = ^v ,2 = = 1 / ■fev ; 

(e) C = diag(l jUiUf); (f) 4 = ^^3 
and ~ 2 3 ! 2 fc ■ Thus, within the uniaxial stress 

framework, the deformation state of the smooth- 
muscle material is defined by the following four 
scalar parameters: X \ , , K and 2 /c . 


2) Stress State 


In the case of the uniaxial tensile test in the e 3 
direction, only the (3,3)-component of different 
stress measures are nonzero. Specifically, the (3,3)- 
component of the first Piola-Kirchhoff stress is 
defined as: 


: h f ( n 3 +n j 


Ms 




kK 

i 


-K(J-1)^ 


A 


3V-^ 

V ^3 J 


2 Phe I ; 




(29) 


where the four terms on the right-hand side of Eq. 
(29) denote respectively the active, the volumetric, 
the steady-state and the viscous stress components, 
I\ = 2/t[ + /t 3 denotes the first invariant (i.e. the 
trace) of ^ . The lateral first Piola-Kirchhoff stresses 
P\ i = P ~i 2 are zero and are defined by the following 
functional relationship: 


P = P 

-m i -*2: 


= k(J- 1)^+^- 


-i/3 r 

3 


PPef 3 


V 

2 h 


3 A, 


A,A V 


(30) 


= 0 


Eq. (30) provides a functional relationship between 
T 3 and f (i.e. f is not an independent 
deformation variable, but rather a variable whose 
value is constrained through the zero-lateral-stress 
condition, Eq. (30)). The corresponding second 
Piola-Kirchhoff stress and Cauchy stress can be 
obtained by combining Eq. (7) with Eqs. (6) and (8). 


3) Evolution Equations 

To complete the uniaxial-stress formulation of the 
smooth-muscle-tissue material model, previously 
defined evolution equations must be used either in 
their original form or after reformulation. 
Specifically: (a) Eq. (1) must be used (in the original 
form) to track the evolution of « 3 and ; (b) Eq. (6) 
must be used (in the original form) to determine 


53 



www.seipub.org/rbb 


Review of Bioinformatics and Biometrics (RBB) Volume 2 Issue 3, September 2013 


temporal evolution of &fc ; and (c) Eq. (28) must be 
used to compute temporal evolution of <^,33 , and 
in turn ^v,33 and cjv,33. 

4) Integration of the Material State: 

To determine the evolution of the material state 
during uniaxial tensile loading, the algebraic and 
differential evolution equations must be combined 
and solved incrementally in time. Loading could be 
either stress-controlled (in which case one would 
determine the temporal evolution of T3 , as well as 
the evolution of the remaining deformation and 
state variables), or stretch-controlled (in which case 
one would compute the temporal evolution of the 
longitudinal stress, as well as the evolution of the 
remaining deformation and state variables). 

Parameterization of the Model 

To parameterize the model overviewed in the 
previous section. Kroon (2010) used isometric 
(constant specimen-length) and isotonic (constant total 
longitudinal stress) uniaxial-tensile-test data for 
swine-carotid-artery smooth-muscle tissue as reported 
in the work of Dillon et al. (1981) and Dillon and 
Murphy (1982). 

Parameterization of the Chemical-state Rate Constants 

In Section III, it was shown that the chemical part of 
the model is defined in terms of seven rate constants, 
fa-fa. To determine these constants. Kroon utilized the 
following approach: (a) since rate constants fa, fa and 
fa are found to be fairly similar for different types of 
smooth muscle, typical values reported in the open 
literature (fa = 4fa, fa = 0.11 s -1 , k7= 0.01 s- 1 ) are adopted; 
(b) for rate constants fa and fa, the values (fc= 0.5 s _1 , fa 
= 0.5 s _1 ) obtained by Hai and Murphy (1988) for 
swine-carotid-artery smooth-muscle are used; (c) as far 
as the remaining rate constants, i.e. fa and fa are 
concerned, it was first recognized that since they 
indirectly or directly affect the (isometric steady-state) 
fraction of cross-bridges in states S3 and sr and thus the 
magnitude of the active stress, these rate constants are 
functions of cytosolic [Ca 2+ ]. As [Ca 2+ ] increases, fa and 
fa start to dominate the remaining rate constants and 
ultimately, at very high values of [Ca 2+ ] solely control 
and result in the maximum isometric steady-state 
value of ()Z3+ m). ki and fa, at a given level of [Ca 2+ ], are 
then obtained by taking advantage of the fact that the 
ratio of the corresponding active stress and that at the 
maximum activation level, is equal to the 


corresponding ratio of the (m + m) sum. This 
procedure yielded fa = fa= 0.14 s _1 for [Ca 2+ ] = 1.6 mM 
(Dillon and Murphy, 1982). 

It should be noted that the values of fa and fa are 
affected by the cytosolic [Ca 2+ ]. The values for these 
two rate constants reported above correspond to the 
cytosolic [Ca 2+ ] of 1.6 mM, the concentration, which 
was maintained constant during the isometric and 
isotonic tensile tests (described in the remainder of this 
section). 

Parameterization of the active part of the model 

As mentioned earlier, the active part of the mechanical 
model is characterized by four parameters: A/, ^3, 
k a and 77. These parameters are determined by fitting 
the model predictions to the swine-carotid-artery 
smooth-muscle tissue isometric tensile test active- 
stress vs. time data obtained by Dillon et al. (1981) and 
Dillon and Murphy (1982). In this type of test, the 
length of the pre-stretched specimen (still under 
tension) is kept constant while extracellular [Ca 2+ ] is 
(abruptly) increased. This causes temporal evolution 
of (and ultimate attainment of a steady-state value for) 
smooth-muscle tissue contraction and the associated 
tensile stress (within the active branch of the 
mechanical model). It should be noted that in 
accordance with Eq. (13), temporal evolution of the 
active stress is affected both by the instantaneous 
values of the chemical variables m-m and by the 
instantaneous value of the mechanical variable tyc . By 
proper time integration of Eq. (1), and using the rate 
constant values listed above, temporal evolution of m- 
J74 can be readily determined. As far as the temporal 
evolution of ^fc is concerned, in accordance with Eq. 
(10), it requires the knowledge of the parameters 77 and 
k 3 . Furthermore, in accordance with Eq. (14), 
determination of the active-stress steady-state value 
entails the knowledge of the parameter . It should 
be noted that the parameter *4 only indirectly, via 
Eqs. (11)-(13), affects the temporal evolution of X fc 
and the active stress during isometric tensile testing. 
Consequently, this parameter was determined using 
the concept of the maximum load-bearing capacity of 
the contractile apparatus. This capacity is defined as 
the maximum stress that the contractile apparatus can 
support, during an abrupt increase in the specimen 
length. Application of a still larger elongation to the 
specimen causes mechanical failure of the contractile 
apparatus (i.e. detachment of all myosin heads from 
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the actin filaments and a complete loss of the load- 
carrying capacity). The value of k a is found by 
matching the model predictions to the experimental 
data pertaining to the maximum load-bearing capacity 
data as obtained in the work of Dillon and Murphy 
(1982). As mentioned above, to obtain the values of the 
four parameters within the active portion of the 
smooth-muscle-tissue mechanical model, the model 
predictions pertaining to the temporal evolution of the 
first Piola-Kirchhoff stress during the activation 
portion of the isometric tensile test, and pertaining to 
the maximum load-carrying capacity obtained in the 
aforementioned test associated with an abrupt sample 
stretch, are matched with their experimental 

counterparts of Dillon et al. (1981) and Dillon and 
Murphy (1982) for isometric testing of swine-carotid- 
artery smooth-muscle tissue. This parameter- 

identification procedure yielded: (a) Mf = 1-91 MPa • (b) 
=1.91 MPa ; (c) /C4 = 0.32 MPa ; and (d) rj = 38.2 
MPa.s. 

Parameterization of the Passive Part of the Model 

As mentioned earlier, the passive part of the smooth- 
muscle-tissue mechanical model is characterized by six 
parameters: K, // !: . ft, a 2 , «?, and r. These parameters 
were determined by fitting model predictions to the 
stretch-rate vs. time data from the isotonic tensile tests 
conducted by Dillon et al. (1981) and Dillon and 
Murphy (1982) on smooth-muscle-tissue from swine 
carotid artery. In this type of test, the specimen is first 
subjected to an isometric test and then the longitudinal 
stress is abruptly decreased (and kept constant 
afterwards). Since the passive part of the mechanical 
model contains elastic elements, the abrupt decrease in 
the stress is associated with an abrupt decrease in the 
corresponding stretch (the so-called " elastic recoil"). 
However, since the passive part of the model also 
includes viscous behavior, the stretch continues to 
evolve with time towards its final (steady-state) value. 
It is the initial abrupt decrease in the stretch, the rate 
of subsequent stretch decrease and the stretch steady- 
state value that are used to determine the 
aforementioned six parameters (via a multiple 
regression curve-fitting procedure). It should be noted 
that within this fitting procedure, the passive portion 
of the smooth-muscle cell tissue is assumed to be 
nearly incompressible, and, hence, K is set to a large 
value (= 10 //;;). The parameter-identification procedure 
described above yielded the following values for the 
six parameters of the passive part of the material 


mechanical model: (a) K = 630 kPa; (b) ju e = 63 kPa; (c) 
P~ 0.14; (d) 02 = 0; (e) 03= 10; and (f) r= 0.50 s. 

Model Validation and Application 

The model presented in the previous section is used in 
this section to reveal the basic response of the swine- 
carotid-artery smooth-muscle tissue during isometric 
and isotonic tensile testing. In addition, model 
predictions are compared with experimental results 
obtained in the work of Dillon et al. (1981) and Dillon 
and Murphy (1982). It should be noted that all the 
calculations carried out in this portion of the 
work involved MATLAB (2012), a general-purpose 
mathematical package. 

Isometric Tensile Test 

The isometric tensile test is conducted using the 
following steps: (a) the test specimen is first stretched 
to the so-called " optimal length" of the smooth-muscle- 
tissue (corresponding to /l 3 =1.6, in the case of swine- 
carotid-artery smooth-muscle tissue (Kamm et al., 
1989)); (b) the pre-stretched specimen is allowed to 
fully relax while, due to a low [Ca 2+ ], the contractile 
apparatus is assumed to be inactive. These 
conditions result in: =1.0, «2 = « 3 = «4 = 0.0, ^fc = ^3 

and A v = J ; (c) the specimen length is then fixed 
throughout the remainder of the test, while 
extracellular [Ca 2+ ] is abruptly increased to a value of 
1.6 mM, Figure 6(a). The [Ca 2+ ] increase promotes the 
bonding of the myosin heads to the actin filaments 
and gives rise to the contraction of the contractile 
apparatus. This in turn causes the development of an 
"active stress and (d) the temporal evolution of the 
(longitudinal component of the) first Piola-Kirchhoff 
stress is then monitored. It should be noted that the 
longitudinal component of the first Piola-Kirchhoff 
stress is numerically equal to the ratio of the 
corresponding force (in the deformed configuration) 
and the specimen cross-sectional area (in the reference 
configuration). 

Temporal evolution of the four chemical state 
variables following an abrupt increase in the 
extracellular [Ca 2+ ] is depicted in Figure 6(b). 
Examination of the results shown in this figure reveals 
that: (a) it takes about 100 s for the chemical-state 
variables to attain the (new) steady-state values; and 
(b) in the new steady-state, most of the cross-bridges 
are in the S4 (latch) state in which they contribute to 
the (active) force transmission. 
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FIGURE 6 EXAMPLES OF THE TYPICAL RESULTS OBTAINED IN THE ISOMETRIC UNIAXIAL TENSILE TEST ON SWINE-CAROTID- 
ARTERY SMOOTH-MUSCLE TISSUE, INVOLVING PRE-STRETCH AND A SUBSEQUENT INCREASE IN THE EXTRA-CELLULAR [Ca 2+ ] 
FROM 0.1 mM TO 1.6 mM. THE RESULTS SHOW TEMPORAL EVOLUTION OF: (a) [Ca 2+ ]; (b) FOUR CHEMICAL-STATE VARIABLES; (c) 
ACTIVE FIRST PIOLA-KIRCHHOFF STRESS; (d) ACTIVE-BRANCH STATUS; (e) MAXIMUM-TO-CURRENT ACTIVE-STRESS RATIO; AND 

(f) FILAMENT CONTRACTION 


Temporal evolution of the longitudinal component of 
the first Piola-Kirchhoff active stress, as predicted by 
the model, is depicted in Figure 6(c). Also shown in 
Figure 6(c) are the experimental results obtained by 
Dillon and Murphy (1982). Examination of the results 
displayed in Figure 6(c) reveals that: (a) the active 


stress also attains its steady-state value after 
approximately 100 s. This finding should be expected 
since the increase in the active stress depends on the 
increase of the four chemical-state variables, n \ ~ n 4 (as 
well as on the evolution of the mechanical state 
variable, -fyc ); and (b) the agreement between the 
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model predictions and the experimental data is 
reasonably good. 

Temporal evolution of the kinematic state of the 
filament within the active branch, as predicted by the 
model, is depicted in Figure 6(d). Examination of this 
figure shows that the contractile apparatus resides, 
during the present isometric tensile test, within the 
contracting state of the filament. Within this state, the 
magnitude of the active stress exceeds the elastic 
resistance, which results in filament contraction, i.e. 

fyc < 0 . Furthermore, the structural integrity of the 
contractile unit is never compromised. 

Temporal evolution of the ratio of the maximum load- 
bearing capacity of the contractile apparatus and the 
current active stress, as predicted by the model, is 
depicted in Figure 6(e). Also shown in Figure 6(e) are 
the experimental results obtained by Dillon and 
Murphy (1982). Examination of the two sets of results 
shows that, with the exception of the short-term 
behavior, the agreement is fairly good. It should be 
noted that for the isometric tensile test, the following 
relationship holds: 

C^act l dX fc = = -(V^) P act,33 

Consequently, the maximum load-carrying capacity 
previously defined by the condition 

~ dY act /dAf c = + K 4 n 4 

can be defined as ^act,max [Xf c / T 3 ) (/c 3 rt 3 +K 4 n 4 )_ 

Temporal evolution of the fiber-contraction stretch 
ratio ^fc , as predicted by the model, is depicted in 
Figure 6(f). Examination of the results shown in this 
figure reveals that it takes about 100 s for this 
mechanical quantity to attain its new steady-state 
value. This finding is expected since the four chemical- 
state variables, active stress and ^f c are mutually 
related, and so are their histories. 


Isotonic Tensile Test 

The isotonic test follows the isometric test and 
includes an abrupt reduction in the longitudinal 
component of the first Piola-Kirchhoff total stress 
(from its value at the end of the isometric test). This 
component of the stress is subsequently held constant 
over the duration of the test. In the present case, the 
longitudinal component of the first Piola-Kirchhoff 
stress is reduced from its initial value of 353.2 kPa to a 
value of 260 kPa, Figure 7(a). 

The temporal evolution of ^3 and A , as predicted by 
the model for the aforementioned abrupt reduction in 
the longitudinal component of the first Piola-Kirchhoff 
stress, is depicted in Figure 7(b). Examination of the 
results depicted in this figure shows that: (a) ^3 drops 
instantaneously (the elastic recoil effect) from the 
initial value of 1.6 to ca. 1.551; (b) 'h continues to 
decrease but at a very low rate; and (c) due to its 
viscous nature, decreases asymptotically towards 
its zero-stress steady-state value. It should be noted 
that the relaxation rate of decreases with an 
increase in r and with a decrease in 

The effect of the longitudinal component of the first 
Piola-Kirchhoff stress level, which is being kept 
constant during the isotonic tensile test, on the initial 
stretch rate of the specimen (following the elastic 
recoil) Aq, as predicted by the model, is depicted in 
Figure 7(c). Also shown in Figure 7(c) are the 
experimental results obtained by Dillon et al. (1981). It 
should be noted that the experimental results of Dillon 
et al. (1981) are used to parameterize the model. Hence, 
the level of agreement observed in Figure 7(c) is a 
measure of the overall quality of the model and the 
success of its parameterization. 
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FIGURE 7. EXAMPLES OF THE TYPICAL RESULTS OBTAINED IN THE ISOTONIC UNIAXIAL-TENSILE-TEST ON SWINE-CAROTID- 
ARTERY SMOOTH-MUSCLE TISSUE, INVOLVING PRE-STRETCH A 3 = 1.6 , FOLLOWED BY AN INCREASE IN THE EXTRA-CELLULAR 
[Ca 2+ ] FROM 0.1 mM TO 1.6 mM, AND CONCLUDED BY A REDUCTION IN THE LONGITUDINAL TENSILE/HOLDING STRESS. THE 

RESULTS SHOW TEMPORAL EVOLUTION OF: (a) THE FIRST PIOLA-KIRCHHOFF P 33 ld STRESS; (b) TEMPORAL EVOLUTION OF 
TOTAL, /b, AND VISCOUS, 2v, LONGITUDINAL STRETCHES; AND THE EFFECT OF THE HOLDING STRESS ON: (c) INITIAL 
REDUCTION-RATE OF /b; AND (d) INSTANTANEOUS VALUE OF te FOLLOWING ELASTIC RECOIL 


The results depicted in Figures 7(a)-(c) all correspond 
to the case of the longitudinal component of the first 
Piola-Kirchhoff stress level being reduced from 353.2 
kPa to 260 kPa. Several other isotonic tests are 
simulated within which the stress was reduced from 
its initial value to different levels. An example of the 
results obtained in these simulations is depicted in 
Figure 7(d). In this figure, the effect of the longitudinal 
component of the first Piola-Kirchhoff isotonic-test- 
stress on the value of obtained after elastic recoil is 
shown. The results displayed in this figure show that, 
as expected, when the stress is kept at its initial value, 
no elastic recoil takes place (^3 remains equal to 1.6). 
On the other hand, as the isotonic stress level 
decreases, the instantaneous drop in O increases (i.e. 
the extent of elastic recoil increases). 

Finite-Element Implementation/Usage of the 
Model 

Finite-Element Implementation 

The material model presented, parameterized and 
applied (to the isometric and isotonic uniaxial-tensile- 
loading test cases) in the previous section is 
implemented, in the present section, into a VUMAT 
user material subroutine and linked with the 
commercial finite-element code ABAQUS/Explicit 
(2011). The task of the user subroutine is to use the 
deformation gradient, current mass density and values 
of the state variables (from the previous time step) in 
order to update the (Cauchy) stress state and the 
material variables at the given material point in the 


current configuration. In the case of the present model, 
the following quantities are defined as state variables: 
m, 112 , m, ru, C v and Afc. 

To calculate/update the stress state and the state 
variables within the new time-step the following 
iterative procedure is used: 

(a) Eq. (1) is integrated in order to update the four 
chemical-state variables m-m; 

(b) The deformation gradients of the three branches 
are set equal to the value of F that is provided to the 
user material subroutine by the ABAQUS/Explicit 
solver; 

(c) Next, the corresponding Jacobian, isochoric part of 
F (Eq. (15)) and Cauchy-Green deformation tensor 
(Eq. (4)) are calculated; 

(d) Since the third branch is purely elastic and, hence, 
its response is time-invariant and dependent on C , 
the volumetric (Eq. (24)) and steady-state (Eq. (25)) 
components of the stress associated with this branch 
are calculated; 

(e) Evolution equation for Aic (Eq. (10)) is then 
integrated over the duration of a given time step in 
order to update this quantity; 

(f) The updated Ate is then used to obtain the active 
part of the stress, i.e. the stress associated with the 
first/active branch (Eq. (14)); 

(e) The evolution equation for the viscous component 
of the stress (Eq. (28)) is combined with the following 
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approximate expression for the time- derivative of S v 


dSy 

dt 


Sy(c /+1 ,c v>I - +1 )-s v (c i+1 ,c v>I .) 


C=const 


At 


(31) 


where subscripts C ;+1 , C v;+1 and C vi are used to 
denote the values of S v computed via Eq. (26) for c 
and C v evaluated at either the current (z+1) or the 
previous ( i ) time increments. The resulting equation is 
solved for the unknown S v (c, +1 ,C iy - +1 ) ; and 


(f) Lastly, the updated value of the viscous stress is 
used to obtain the corresponding updated value of C v , 
i.e. C vJ+l via Eq. (26). It should be noted that Eq. (26) 
represents a system of six linear algebraic equations 
with six unknowns. 


Implementation Validation 

The model implemented into the user-material 
subroutine is then validated by re-simulating the 
isometric and isotonic uniaxial tensile tests previously 
described in Section VII. The results obtained are 
identical to their counterparts reported in Figures 6(a)- 
(f) and 7(a)-(d), confirming the validity of the model 
implementation. 

Application of the Model to a Simple Contraction Problem 

1) Problem Formulation 

In this section, the user-material subroutine 
described and validated in the previous two 
sections is utilized in the context of a simple 
smooth-muscle contraction problem in order to 
demonstrate its utility. The problem analyzed 
involves a tubular structure made of smooth 
muscles with the cells' longitudinal direction 
aligned with the circumferential direction of the 
structure. In addition, tube intra-luminal content is 
modeled as a prolate spheroidal with the major/ 
unique axis collinear with the tube axis. The minor 
axis of the spheroidal, on the other hand, is made 
slightly larger than the inner radius of the tube (i.e. 
the spheroidal and the tube are initially in contact 
over a small area). The contact stresses activate 
smooth muscles, in the contact region, by 
increasing the cytosolic [Ca 2+ ]. This, in turn, causes 
longitudinal contraction of the activated 
("awakened'') smooth-muscle cells, which gives rise 
to the local circumferential/radial contraction of the 
tube. The contraction of the tubular structure then 
increases the ellipsoidal/tube contact stresses and, 
hence, the net force acting on the ellipsoidal in the 
tube-axis direction. This results in the forward 
displacement of the ellipsoidal. As a result of the 


downward displacement of the ellipsoidal, a new 
spheroidal/tube contact region is generated and a 
new tubular section of the smooth muscle is 
activated. This, in turn, results in further forward 
advancement of the ellipsoidal, and the whole 
sequence of processes described above repeats. 

2) Geometrical Model 

The initial geometry of the computational domain, 
consisting of the previously described tubular and 
ellipsoidal structures, is depicted in Figure 8(a). 
The initial dimensions of the two structures are as 
follows: (a) tube - length = 30 mm; outer radius = 5 
mm; inner radius = 3 mm; and (b) ellipsoid - major 
axis = 6 mm; minor axis = 3.025 mm. The 
longitudinal axes of the tube/ellipsoidal are aligned 
in the z-direction. 



FIGURE 8. TYPICAL: (a) GEOMETRICAL; AND (b) MESHED 
MODELS OF THE TUBE AND ELLIPSOIDAL USED IN THE 
SIMULATIONS OF SMOOTH-MUSCLE CONTRACTION 

3) Meshed Model 

The geometrical model described above is meshed 
in the following way: (a) the tube is meshed using 4, 
40 and 60 elements in the through-the-thickness, 
circumferential and longitudinal directions, 
respectively. The elements used are eight-node, 
hexahedral, first-order, continuum elements; and (b) 
the ellipsoid is meshed using 616 quadrilateral and 
56 triangular membrane elements. For simplicity, 
the ellipsoidal is treated as a rigid body with six 
(three translational and three rotational) degrees of 
freedom. Typical finite-element meshes used in 
the present work for the two parts of the 
computational domain are displayed in Figure 8(b). 
All meshing carried out in the present work was 
done using HyperMesh (2012), a general-purpose 


59 





www.seipub.org/rbb 


Review of Bioinformatics and Biometrics (RBB) Volume 2 Issue 3, September 2013 


pre-processing software from Altair Engineering, 
Inc. 

4) Material Model 

The tubular structure is assigned the smooth- 
muscle material model overviewed, parameterized 
and validated in the previous sections. As far as 
the spheroidal-shaped intra-luminal content is 
concerned, it is treated as a rigid body and, hence, 
no material model formulation (except for the mass 
density) had to be specified. 

5) Initial Conditions 

The tubular structure, in its original configuration, 
is assumed to be in equilibrium with the interacting 
ellipsoidal and to contain stresses developed 
within the passive portion of the material as a 
result of the tube/ellipsoidal contact stresses. Both 
the tube and the ellipsoidal are initially assumed to 
be quiescent. 

6) Boundary Conditions 

Fixed-displacement boundary conditions are applied 
in the longitudinal direction of the tube, to one of 
the tube ends. In addition, zero-stress conditions 
are applied to all the free surfaces of the tube and 
the ellipsoidal. 

7) Contact Interactions 

The tube-ellipsoidal normal interactions are 
analyzed using a penalty-contact algorithm. Within 
this algorithm, (normal) penetration of the 
contacting surfaces is resisted by a set of linear 
springs which produce a contact pressure that is 
proportional to the depth of tube/ellipsoidal 
penetration. To minimize penetration, the maximum 
values, which still ensure computational stability, 
are assigned to the (penalty) spring constants. 
Force equilibrium in a direction collinear with the 
contact-interface normal then causes the 
penetration to acquire an equilibrium (contact- 
pressure-dependent) value. It should be noted that 
no contact pressures are developed unless (and 
until) the nodes on the (ellipsoidal) "slave surfaces" 
contact/penetrate the (tube inner) "master surfaces". 
On the other hand, the magnitude of the contact 
pressure that can be developed is unlimited. As far 
as the tangential tube-ellipsoidal interactions 
(responsible for transmission of the shear stresses 
across the contact interface) are concerned, they are 
modeled using the Coulomb friction law. Within 
this law, the maximum value of the shear stresses 
that can be transmitted (before the contacting 


surfaces begin to slide) is defined by a product of 
the contact pressure and a static (before sliding) 
and a kinetic (during sliding) friction coefficient. 

8) Computational Procedure 

Explicit, transient, non-linear-dynamics finite 
element analyses are employed while ensuring that 
the stability criterion is met through the proper and 
adaptive selection of the time increment. As 
mentioned earlier, all the calculations carried out in 
the present work were performed using ABAQUS/ 
Explicit (2011). 

9) Computational Accuracy, Stability and Cost 

A standard mesh sensitivity analysis was carried 
out (the results not shown for brevity) in order to 
ensure that the results obtained are accurate, i.e. 
insensitive to further reductions in the size of the 
elements used. A typical 10-second computational 
analysis followed by a detailed post-processing 
data reduction analysis required on average 30 
minutes of (wall-clock) time on a 12-core, 3.0 GHz 
machine with 16 GB of memory. 

10) Typical Results 

An example of the typical results is depicted in 
Figures 9(a)-(d). The results depicted in these 
figures show the spatial distribution of the two 
materials/bodies, i.e. of the smooth-muscle tissue in 
the form of a tubular section, and the ellipsoidal- 
shaped intra-luminal content at four times (2 sec, 4 
sec, 6 sec, and 8 sec) measured from the instant of 
initial activation of the muscle. Examination of the 
results shown in Figures 9(a)-(d) clearly reveals 
that: (a) the contact stresses at the tube/ellipsoidal 
interfaces promote smooth-muscle tissue contraction; 
(b) as a result of this contraction, forward 
momentum is imparted to the ellipsoidal, giving 
rise to its advancement along the tube axis; and (c) 
once the ellipsoidal has advanced and the 
contracted smooth-muscle tissue no longer 
experiences contact pressure, the tissue relaxes and 
the tube acquires its original cross-section. 



60 



Review of Bioinformatics and Biometrics (RBB) Volume 2 Issue 3, September 2013 


www. seipub . org/rbb 



FIGURE 9. SPATIAL DISTRIBUTION OF THE SMOOTH-MUSCLE 
TISSUE IN THE FORM OF A TUBULAR SECTION, AND THE 
ELLIPSOIDAL-SHAPED INTRA-LUMINAL CONTENT, AT FOUR 
TIMES MEASURED FROM THE INSTANT OF INITIAL 
ACTIVATION OF THE MUSCLE: (a) 2 sec; (b) 4 sec; (c) 6 sec; AND 
(d) 8 sec 

Temporal evolution of the axial position of the 
ellipsoidal center of mass is shown in Figure 10(a). 
It is seen that the motion of the ellipsoidal is not 
smooth but rather jerky. That is, occasionally the 
ellipsoidal "runs away" from the contracting 
smooth-muscle tissue and comes to rest. Then, the 
contact stresses activate the tissue in the region 
surrounding the contact interface, causing additional 
advancement of the ellipsoidal. Temporal evolution 
of the radial displacement at one of the points 
which is located at the interior surface of the tube, 
and halfway along its length, is depicted in Figure 
10(b). It is seen that as expected: (a) the smooth- 
muscle tissue is initially in the relaxed state, which 
is characterized by zero radial displacement at the 
location described above; (b) as the ellipsoidal 
begins to approach the target point, the radial 
displacement first begins to acquire positive values 
(i.e. the tube begins to expand in the radial 
direchon in order to accommodate the slightly- 
oversized ellipsoidal); (c) as the tube/ellipsoidal 
contact interface passes through the target point. 


the radial displacement begins to decrease and 
ultimately acquires negative values (as a result of 
the smooth-muscle tissue contraction); and (d) once 
the ellipsoidal has advanced sufficiently far from 
the target point, the muscle-tissue begins to relax, 
causing the radial displacement to start approaching 
its initial (zero) value. 




FIGURE 10. (a) TEMPORAL EVOLUTION OF THE AXIAL 
POSITION OF THE ELLIPSOIDAL CENTER OF MASS; AND (b) 
TEMPORAL EVOLUTION OF THE RADIAL DISPLACEMENT AT 
ONE OF THE POINTS WHICH IS LOCATED AT THE INTERIOR 
SURFACE OF THE TUBE, AND HALFWAY ALONG ITS LENGTH 

Summary and Conclusions 

Based on the results obtained in the present work, the 
following summary remarks and main conclusions can 
be drawn: 

1. A critical assessment, parameterization (using 
experimental data, for isometric (constant 
length) and isotonic (constant holding stress) 
uniaxial-tensile-tests of swine carotid artery. 
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obtained by Dillon et al. (1981) and Dillon and 
Murphy (1982)) and validation of the smooth- 
muscle-tissue material model originally 
proposed by Kroon (2010) has been conducted. 
It has been shown that, for the most part, the 
model correctly accounts for key aspects of the 
smooth-muscle-tissue behavior under these 
loading conditions. 

2. The model is next implemented into a user- 
material subroutine within which, for a given 
incremental deformation, the stress state as 
well as the values of the material-state 
variables are updated over a small time step. 

3. The user-material subroutine is then linked 
with a commercial finite-element program and 
used to carry out a computational analysis of a 
simple smooth-muscle contraction problem 
involving: (a) a tubular structure surrounded 
with smooth-muscle tissue (with muscle-cell 
longitudinal directions aligned in the tube 
circumferential direction); and (b) an intra- 
luminal content in the form of a prolate 
spheroidal with the spheroidal major axis 
coinciding with the tube axis. 

The results of the computational investigation have 
clearly demonstrated how activation/ contraction of the 
smooth-muscle tissue triggered by the tube/spheroidal 
contact stresses results in a forward momentum/ 
motion of the intra-luminal content. 
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